#include <iostream>
#include <fstream>
#include <sstream>
#include <ctime>
#include <ratio>
#include <chrono>
#include <omp.h>
#include "mpi.h"

using namespace std;

const int col = 8399;      // 列数
const int elinenum = 4535; // 被消元行数
const int num_thread = 4;  // OpenMP线程数
const int mpisize = 8;     // MPI进程数
const int bytenum = (col - 1) / 32 + 1; // 每个实例中的byte型数组数

bool parallel = 1;
int isupgrade;
int tmp = 0;

class bitmatrix {
public:
    int mycol;    // 首项
    int *mybyte;
    bitmatrix() { // 初始化
        mycol = -1;
        mybyte = new int[bytenum];
        for (int i = 0; i < bytenum; i++)
            mybyte[i] = 0;
    }
    bool isnull() { // 判断当前行是否为空行
        return mycol == -1;
    }
    void insert(int x) { // 数据读入
        if (mycol == -1) mycol = x;
        int a = x / 32, b = x % 32;
        mybyte[a] |= (1 << b);
    }
    void doxor(bitmatrix b) { // 两行做异或操作，由于结果留在本实例中，只有被消元行能执行这一操作，且异或操作后要更新首项
        for (int i = 0; i < bytenum; i++)
            mybyte[i] ^= b.mybyte[i];
        for (int i = bytenum - 1; i >= 0; i--)
            for (int j = 31; j >= 0; j--)
                if ((mybyte[i] & (1 << j)) != 0) {
                    mycol = i * 32 + j;
                    return;
                }
        mycol = -1;
    }
};

bitmatrix *eliminer = new bitmatrix[col], *eline = new bitmatrix[elinenum], *pass = new bitmatrix[mpisize];

void readdata() {
    ifstream ifs1("/mpi/test/1.txt");
    string temp;
    while (getline(ifs1, temp)) {
        istringstream ss(temp);
        int x;
        int trow = 0;
        while (ss >> x) {
            if (!trow) trow = x; // 第一个读入元素代表行号
            eliminer[trow].insert(x);
        }
    }
    ifs1.close();

    ifstream ifs2("/mpi/test/2.txt");
    int trow = 0;
    while (getline(ifs2, temp)) {
        istringstream ss(temp);
        int x;
        while (ss >> x) {
            eline[trow].insert(x);
        }
        trow++;
    }
    ifs2.close();
}

void printres() { // 打印结果
    for (int i = 0; i < elinenum; i++) {
        if (eline[i].isnull()) {
            puts("");
            continue;
        }
        for (int j = bytenum - 1; j >= 0; j--) {
            for (int k = 31; k >= 0; k--) {
                if ((eline[i].mybyte[j] & (1 << k)) != 0) {
                    printf("%d ", j * 32 + k);
                }
            }
        }
        puts("");
    }
}

void dowork() { // 消元
    int rank;
    MPI_Status status;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank); // 获取当前进程号
    int r1 = rank * (elinenum / mpisize);
    int r2 = (rank == mpisize - 1) ? elinenum - 1 : (rank + 1) * (elinenum / mpisize) - 1;
    double st = MPI_Wtime(); // 计时开始

    #pragma omp parallel if(parallel), num_threads(num_thread), private(rank)
    {
        for (int i = col - 1; i >= 0; i--) {
            if (!eliminer[i].isnull()) {
                #pragma omp for
                for (int j = r1; j <= r2; j++) {
                    if (eline[j].mycol == i) {
                        eline[j].doxor(eliminer[i]);
                    }
                }
            } else {
                isupgrade = -1;
                int t = -1;
                if (rank != 0) {
                    #pragma omp single
                    {
                        for (int k = r1; k <= r2; k++) {
                            if (eline[k].mycol == i) {
                                pass[rank] = eline[k];
                                t = k;
                                MPI_Send(&t, 1, MPI_INT, 0, k + elinenum * 4 + 3, MPI_COMM_WORLD);
                                MPI_Send(&pass[rank].mybyte[0], bytenum, MPI_INT, 0, k + 3, MPI_COMM_WORLD);
                                MPI_Send(&pass[rank].mycol, 1, MPI_INT, 0, k + elinenum + 3, MPI_COMM_WORLD);
                                break;
                            }
                        }
                        if (t == -1) {
                            MPI_Send(&t, 1, MPI_INT, 0, k + elinenum * 4 + 3, MPI_COMM_WORLD);
                            MPI_Send(&pass[rank].mybyte[0], bytenum, MPI_INT, 0, k + 3, MPI_COMM_WORLD);
                            MPI_Send(&pass[rank].mycol, 1, MPI_INT, 0, k + elinenum + 3, MPI_COMM_WORLD);
                        }
                    }
                } else {
                    #pragma omp single
                    {
                        for (int k = mpisize - 1; k > 0; k--) {
                            MPI_Recv(&t, 1, MPI_INT, k, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                            MPI_Recv(&pass[k].mybyte[0], bytenum, MPI_INT, k, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                            MPI_Recv(&pass[k].mycol, 1, MPI_INT, k, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                            if (t != -1) {
                                eliminer[i] = pass[k];
                                isupgrade = t;
                            }
                        }
                        for (int k = r1; k <= r2; k++) {
                            if (eline[k].mycol == i) {
                                eliminer[i] = eline[k];
                                isupgrade = k;
                                break;
                            }
                        }
                        for (int k = 1; k < mpisize; k++) {
                            int t1 = k * (elinenum / mpisize), t2 = (k == mpisize - 1) ? elinenum - 1 : (k + 1) * (elinenum / mpisize) - 1;
                            MPI_Send(&isupgrade, 1, MPI_INT, k, 0, MPI_COMM_WORLD);
                            if (isupgrade != -1 && t2 >= isupgrade) {
                                MPI_Send(&eliminer[i].mybyte[0], bytenum, MPI_INT, k, 1, MPI_COMM_WORLD);
                                MPI_Send(&eliminer[i].mycol, 1, MPI_INT, k, 2, MPI_COMM_WORLD);
                            }
                        }
                    }
                }
                if (rank != 0) {
                    MPI_Recv(&isupgrade, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
                    if (isupgrade != -1 && r2 >= isupgrade) {
                        MPI_Recv(&eliminer[i].mybyte[0], bytenum, MPI_INT, 0, 1, MPI_COMM_WORLD, &status);
                        MPI_Recv(&eliminer[i].mycol, bytenum, MPI_INT, 0, 2, MPI_COMM_WORLD, &status);
                    }
                }
                if (isupgrade != -1 && r2 >= isupgrade) {
                    #pragma omp for
                    for (int j = r1; j <= r2; j++) {
                        if (eline[j].mycol == i && j != isupgrade) {
                            eline[j].doxor(eliminer[i]);
                        }
                    }
                }
            }
        }
    }

    if (rank != 0) {
        for (int k = r1; k <= r2; k++) {
            MPI_Send(&eline[k].mybyte[0], bytenum, MPI_INT, 0, k + 3 + elinenum * 2, MPI_COMM_WORLD);
            MPI_Send(&eline[k].mycol, 1, MPI_INT, 0, k + 3 + elinenum * 3, MPI_COMM_WORLD);
        }
    } else {
        for (int k = 1; k < mpisize; k++) {
            int t1 = k * (elinenum / mpisize), t2 = (k == mpisize - 1) ? elinenum - 1 : (k + 1) * (elinenum / mpisize) - 1;
            for (int l = t1; l <= t2; l++) {
                MPI_Recv(&eline[l].mybyte[0], bytenum, MPI_INT, k, l + 3 + elinenum * 2, MPI_COMM_WORLD, &status);
                MPI_Recv(&eline[l].mycol, 1, MPI_INT, k, l + 3 + elinenum * 3, MPI_COMM_WORLD, &status);
            }
        }
        double ed = MPI_Wtime(); // 计时结束
        cout << "parallel: " << ed - st << endl;
    }
}

int main(int argc, char *argv[]) {
    readdata();
    MPI_Init(&argc, &argv);
    dowork();
    MPI_Finalize();
    // printres();
    return 0;
}
